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Abstract: 



We construct dynamical models for a sample of 36 nearby galaxies with Hubble Space 
Telescope photometry and ground-based kinematics. The models assume that each galaxy 
is axisymmetric, with a two-integral distribution function, arbitrary inclination angle, a 
position-independent stellar mass-to-light ratio T, and a central massive dark object (MDO) 
of arbitrary mass M.. They provide acceptable fits to 32 of the galaxies for some value of 
M. and T; the four galaxies that cannot be fit have kinematically decoupled cores. The 
mass-to-light ratios inferred for the 32 well-fit galaxies are consistent with the fundamental 
plane correlation T oc L° 2 , where L is galaxy luminosity. In all but six galaxies the models 
require at the 95% confidence level an MDO of mass M. ~ 0.006M bu i ge = 0.006TL. Five of 
the six galaxies consistent with M. = are also consistent with this correlation. The other 
(NGC 7332) has a much stronger upper limit on M. . We consider various parameteriza- 
tions for the probability distribution describing the correlation of the masses of these MDOs 
with other galaxy properties. One of the best models can be summarized thus: a fraction 
/ ~ 0.97 of galaxies have MDOs, whose masses are well described by a Gaussian distribution 
in log(M./M bu i ge ) of mean -2.27 and width ~ 0.07. 



1 Introduction 

The evidence that massive dark objects (MDOs) are present in the centers of nearby galaxies is reviewed by 
Kormendy & Richstone (1995; hereafter KR95). Further evidence that post-dates this review is described 
by Bender, Kormendy & Dchnen (1997), van der Marel et al. (1997) and Kormendy et al. (1997a). The 
MDOs are probably black holes, since star clusters of the required mass and size are difficult to construct 
and maintain, and since black-hole quasar remnants are expected to be common in galaxy centers; however, 
this identification is not important for the purposes of this paper. Following Kormendy (1993a), KR95 
suggest that at least 20% of nearby hot galaxies (ellipticals and spiral bulges) have MDOs and point out 
that the observed MDO masses exhibit the correlation M, ~ 0.003Mb u i go , where Mb u i g c is the mass of the 
hot stellar component of the galaxy. (Throughout this paper we use the word "bulge" to refer to the hot 
stellar component of a galaxy, whether elliptical or spiral.) For a "bulge" with constant mass-to- light ratio 
T and luminosity L, Mb u i g o = TL. 

The machinery for modelling the kinematics of hot galaxies to determine whether MDOs are present has 
increased steadily in sophistication over the past two decades. The earliest models (e.g. Young et al. 1978) 
fitted only the line-of-sight dispersion of spherical galaxies and assumed that the stellar distribution function 
was isotropic. Modern programs (e.g. Rix et al. 1997; Gebhardt et al. 1997) fit the entire line-of-sight velocity 
distribution for arbitrary axisymmetric galaxy models. While the most general and accurate possible models, 
and the highest resolution spectroscopic observations, were needed to establish the presence of the first few 
MDOs, we have learned with experience that estimates of the MDO mass based on cruder models and 
observations are usually fairly accurate. An example is the MDO in M87: Young et al. (1978) estimated 
the mass to be ~ 5 x 10 9 M Q from spherical, isotropic models, very close to the 3 x 10 9 M Q determined by 
Harms et al. (1994) from HST spectra of a ring of ionized gas at 20 pc from the center. 

This experience suggests that it is worthwhile to estimate MDO masses using relatively simple models applied 
to a large sample of galaxies. We cannot yet insist on HST spectroscopy for our sample, since this is still 
available only for a few galaxies; on the other hand HST photometry is available for over 60 hot galaxies. In 
this paper we examine a sample of 36 hot galaxies for which both HST photometry and reasonable quality, 
ground-based, long-slit spectroscopy are available. We look for evidence of MDOs among these by fitting 
two-integral axisymmetric dynamical models to the data for each galaxy. These are not the most general 
types of models, but they are quick to compute and will guide us towards galaxies to which we should 
apply more precise (and expensive) observations and models. Our results also provide a first look at the 
statistical distribution of MDOs as a function of galaxy luminosity and other parameters. They do not 
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establish unambiguously that an MDO is present in any individual galaxy. 

The paper is organized as follows. The next section gives a brief outline of the data we use. This is followed 
by a detailed description of our modelling procedure and the assumptions that go in to it. Section 4 presents 
results for individual galaxies. What these tell us about the MDO mass distribution is tackled in Section 5. 
Finally Section 6 sums up. 



2 Data 

Our sample consists of all reasonably dust-free hot galaxies with HST photometry and ground-based velocity 
dispersion and rotation velocity profiles. The sample contains 36 galaxies, listed in Table 1 along with details 
of the sources of the observations we use. The Appendix contains comments about some of the galaxies. 
Most of the objects were observed prior to the first HST servicing mission with the Planetary Camera (0.043" 
per pixel) through filter F555W (roughly Johnson V). The sample includes galaxies observed in a number 
of HST programs (Lauer et al. 1992a, b, Grillmair et al. 1994, Jaffe et al. 1994, Forbes et al. 1995, Lauer et 
al. 1995); the reduction procedures are described by Lauer et al. (1995), Byun et al. (1996), and Faber et 
al. (1997). Some of the galaxies show evidence for nuclear activity. For each of these the table lists a radius 
fimin, inside which non-stellar radiation probably makes a significant contribution to the observed light. 

The HST data extend only to about 10" from the centre. For most of the galaxies we take published 
ground-based photometry and join it smoothly to the HST data to obtain a global photometric profile. For 
the remaining galaxies, we assume that the outer parts are well described by an i? 1 / 4 profile with the same 
flattening as the outermost HST isophote. We take the effective radius of the i? 1 / 4 profile from the literature, 
if available; otherwise we estimate it by fitting to the HST photometry. 

Table 1 also lists the sources for our kinematical data. We restrict ourselves to reasonable quality CCD-based 
spectroscopy and do not use kinematical data beyond about two-thirds of the maximum radius for which 
photometry is available. Whenever there are many sources for a given slit position of a galaxy, we generally 
choose those with the best seeing. If an estimate for the seeing is unavailable we simply assume a FWHM 
of 2" : the MDO masses yielded by our models are fairly insensitive to the precise value used, as long as it 
lies between 1" and 3" . 

The observations yield line-of-sight rotation speeds and velocity dispersions, convolved with seeing and 
averaged over spatial "bins" determined by the slit width and pixel size. We combine the measured rotation 
speed Vj and the velocity dispersion <jj in each bin j to obtain an estimate of the second-order moment 
$ = v j a j ' wmcn i s tne i n P u t used by our models. Strictly speaking, the Vj and cfj quoted by observers 
do not individually have any direct connection with the moments of the line-of-sight velocity profiles (VPs) , 
since they are usually obtained by fitting Gaussians to the VPs (van der Marel & Franx 1993). However, 
tests with flattened isotropic toy galaxies (Dehnen & Gerhard 1994; Magorrian & Binney 1994) show that 
there is typically an almost-constant difference of about 10% between the combination Vj + a? and the true 
second-order moments, with the sign of the difference changing from the major to the minor axis. Figure 1 
shows a typical example. Since these differences are almost constant, they do not affect the MDO masses 
fitted by our models. They could, however, have a small effect on the fitted mass-to-light ratios. 

Deriving the observational uncertainty A/j,j is a vexing problem, because analysis methods used by different 
observers yield a range of error estimates Avj and Aaj, which usually do not take systematic effects, such as 
template mismatch, into account. We have tried the following three methods of dealing with this problem: 

(i) Simply take all quoted errors at face value; 

(ii) Replace Avj with max(A^,5kms _1 ) and similarly for Act-,-; 

(iii) Scale the errors for each exposure along each slit position such that they are consistent with axisymmetry. 
More precisely, suppose there are n measurements cr+ along one side of a galaxy, with corresponding 
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measurements along the other side. We scale the Act, by a constant factor such that 

X >-Y =n (1) 

and similarly for the Avj (Davies & Birkinshaw 1988). 

Notice that the last method implicitly assumes that the errors are Gaussian. This is almost certainly wrong, 
but it is the best we can do given the heterogeneous nature of our data. Given the "improved" observational 

1 /2 

errors, the error in fij is A/ij = (vj(Avj) 2 + ctJ(AcTj) 2 ) / fj,j to first order. We use (iii) above wherever 
possible, but the results of our models are usually not significantly affected by which procedure we employ. 



3 Modelling Procedure 

We assume that each galaxy is axisymmetric with some unknown inclination angle i, and work in cylindrical 
coordinates (R, <p, z) where the z-axis is the symmetry axis of the galaxy. A lower bound on i comes from 
requiring that all isophotcs have an intrinsic axis ratio no less than 0.3 (i.e. no flatter than E7). Each 
galaxy can have a central MDO of arbitrary mass, but otherwise the mass-to-light ratio T is assumed to 
be independent of position. The distribution function of the stars is assumed to be a function only of 
two integrals of motion, the energy and the z-componcnt of angular momentum. The advantage of these 
assumptions is that the (even part of the) kinematics follows uniquely from the three-dimensional luminosity 
distribution v(R, z) (e.g., Lynden-Bell 1962; Dejonghe 1986). The disadvantage is that there is no reason why 
real galaxies should obey our assumptions. In particular, bright, core galaxies are usually non-rotating, and 
many studies of them (e.g., van der Marel 1991) show some evidence for radial anisotropy. Our two-integral 
models of flattened, non-rotating galaxies are tangentially anisotropic. 

We use the modelling procedure introduced by Binney, Davies & Illingworth (1990) to predict the kinematics 
of each galaxy for any assumed inclination angle i and MDO mass M # . It predicts the second-order moments, 
convolved with seeing and averaged over the same j — 1, . . . , n spatial bins used in the observations. The 
procedure is as follows: 

1. Use a scheme based on maximum penalized likelihood to find a smooth luminosity density v{R, z) that 
projects to an acceptable fit to the observed surface brightness (Magorrian 1997). The density v is not 
uniquely determined by the surface brightness unless the galaxy is edge-on (Rybicki 1987). Romanowsky 
& Kochanek (1997) demonstrate that even for quite high inclinations there can be a large range in v 
consistent with a given surface brightness; however, they find that the range of projected second-order 
moments associated with this uncertainty is quite small. We have carried out some experiments that 
confirm that the allowable MDO masses are not strongly affected by the indeterminacy in v(R, z); 

2. Calculate the gravitational potential and forces using an assumed stellar mass-to-light ratio T and 
MDO mass M.; 

3. Use the Jeans equations to calculate the second-order moments vv^ and vv 2 R — vv 2 z \ 

4. Project the luminosity-weighted zeroth- and second-order moments of the line-of-sight velocity along 
the line of sight; convolve with seeing; and average over the same spatial bins used in the observations. 

Dividing the binned, seeing-convolved second-order moment by the corresponding zeroth-order one yields 
the model's predictions /t|(i,T ,M.) in each bin. These predictions scale trivially with mass-to-light ratio 
T through 

T, M.) = -^A 2 (i, T , TM./To) • (2) 

T 
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3.1 Estimation of m. and r. 



We assume that the measurement errors in the fij are Gaussian and uncorrelated. Then the likelihood of 
the photometric and kinematic data D given the model parameters (i, T, M.) is 

pph,T,M.)cxexp(-ix 2 ), (3) 

where 



i=i 



x!(i , T ,M. )s £(E^) 2 . (4) 



We obtain the best-fitting values M. and T and their confidence intervals as follows. By Bayes' theorem, 
the posterior distribution of i, T and M, given the data D is 

p(i, T, M. | £>) oc p(D | T, M.) p(i | </) p(T) p(M.), (5) 

where we have made the assumption that i, T and M. are a priori independent. Our priors p(M,) and p(T) 
are flat in M, and log T respectively. We make the reasonable assumption that the prior for i depends only 
on the observed axis ratio q' of the galaxy. Then p(i | q') can be related to N(q) dq, the probability that a 
randomly chosen galaxy will have an intrinsic axis ratio lying between q and q + dq, by a further application 
of Bayes' theorem: 

p(Qp(g' I i) 

«p(i) J P(q' I i,q)N(q)dq (6) 



P(i I <?') = 



cx TV 



(g^-cos 2 *) 1 ^ 



v/g' 2 - cos 2 i V sin i 

where we have made the natural assumption that p(i) = sini and have used the relation p(q' \ i,q) oc 
q'5(q 2 sin 2 i + cos 2 i — q' 2 ). We approximate the N(q) obtained by Tremblay & Merritt (1995) by a Gaussian 
centred on q = 0.7 with standard deviation 0.1. Our results are only very weakly dependent on this form. 

We are interested mainly in M, and T, not in i. Marginalizing the latter, we get the joint posterior 
distribution of T and M. as 

P (T, M. | D) = v{D | T, M.) p(T) p(M.), (7) 

where 

p(£> | T, M.) =J P (D\ i, T, M.) p(i | </) di. (8) 
The posterior distributions p(M. | _D) and p(T | D) follow by marginalizing (7) again. 
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4 Results for individual galaxies 



We have made models of each of our 36 galaxies for a range of M, > 0, T and i. The models do not 
provide adequate descriptions of the kinematics of four of the galaxies (NGC 1700, NGC 4365, NGC 4494 
and NGC 4589). For these four, Figure 2(a) shows how x 2 of equation (4) varies with MDO mass M, and 
inclination angle i. Figure 2(b) plots the kinematics of the models with the best-fitting values of M, against 
the observations. All four galaxies are known to have kinematically distinct cores (e.g., Forbes et al. 1996), 
so it is perhaps not surprising that our axisymmetric models do not work for them. We omit these four in 
the demographical analysis in the next section. For comparison, only two (NGC 3608 and NGC 4278) of the 
32 galaxies that our models do fit are known to have kinematically distinct cores. 

The models describe the kinematics of all of the remaining 32 galaxies reasonably well for some value of M, . 
Figure 3(a) shows the posterior distribution p(T, M, \ D) in each case. Figure 3(b) shows the kinematics of 
the best-fitting models along each slit position fitted. Six of these galaxies have independent determinations 
of the mass of a central MDO. The comparison between the best-fit masses as determined here and the mass 
estimates in the literature for these six is presented in Figure 4. For all but one galaxy we obtain MDO 
masses that are in good agreement with those from earlier work. This gives us some degree of confidence 
in the assumptions that go into our models. The one exception is NGC 3115 for which Kormendy et al. 
(1996a) claim an MDO mass of about 2 x 1O 9 M , some four times larger than our present mass estimate. 

Our models imply that only three of the 32 galaxies (NGC 2778, NGC 4467, and NGC 7332) are consistent 
(at the 68% confidence level) with M, = 0. However, Figures 3(a) and (b) show that the available data for 
each of these three are also consistent with a reasonably large value of M,. NGC 7332 has the strongest 
upper limit on M.. In fact, the central dip in its dispersion profile is suggestive of either a mass-to-light ratio 
that decreases close to the centre, or else strong tangential anisotropy. Kormendy (1993b) has suggested that 
its formation history may be different from the other galaxies in the sample. We do not, however, omit it 
in the analysis below. There are a further three galaxies (NGC 4168, NGC 4473 and NGC 4636) consistent 
with M, — at the 95% confidence level. All the rest have M, > 0. 

Seven of the 32 galaxies show evidence for nuclear activity or strong dust obscuration. For these we make 
two types of models: one under the naive assumption that all the observed light near the galaxy centre 
comes from stars, the other that only uses the photometry beyond a radius i? m i n , where i? m i n is given in 
Table 1. We find that the MDO masses predicted by the two types of models generally agree quite well. 
This is unsurprising given the relatively poor spatial resolution of the kinematical data. In what follows we 
use only the MDO masses obtained by omitting photometry within i? m i n . 

Table 2 lists the 68% confidence bounds that our models place on M, and T for the 32 galaxies. The 
correlations between T and L and between M, and Mb u i ge are plotted on Figure 5. Ignoring the error bars 
on T the formal best-fit straight line to (log L, log T) is 

log(T fit /T ) = (-1.12 ± 0.33) + (0.18 ± 0.03) Iog(L/L ) (9) 

with an RMS deviation between logT nt and logT of 0.12. Our crude fit is broadly consistent with the 
fundamental-plane correlation T oc L 2 predicted using the virial theorem (e.g., Faber et al. 1987; Bender, 
Burstein & Faber 1992). Similarly, the correlation between M. and Mb u i ge for those galaxies with M, > 
can be described by 

Iog(M., fit /M ) - (-1.82 ± 1.36) + (0.96 ± 0.12) log(M bulgc /M ), (10) 

with an RMS (logM.^t — logM.) of 0.49. This result is consistent with the proportionality M, oc Mbuige 
that was first pointed out by Kormendy (1993a) and KR95. This apparent correlation is the subject of the 
next section. 

Finally, we check whether there is any correlation between the residuals xi = logM.^ — logM. jnt ,i and 
Hi = logTj — log Yet,;. One might expect a negative correlation if our models were fitting spuriously high 
MDO masses for some galaxies, thus depressing the fitted value of T. Figure 5(c) shows that there is no 
such correlation. The correlation coefficient r xy — 0.014, which is not significant. 
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Table 3. The five parametcrizations for p(x \ ui) considered here. The variable x = M./M bulgo where M. is the mass of the 
MDO and M bu igo is the mass of the hot stellar component of the galaxy. For a given set of parameters u>, the probability that 
a galaxy has an MDO with mass in the range [x, x + dx] is p(x \ ui, P) dx = (1 — f)8(x) dx + / p + (x \ d, P) dx, where / is the 
fraction of galaxies with M. > and the N(ui) in p+(x \ ui,P) is a normalizing factor (equation (12)). The prior probability 
p(ui | P) is assumed to be flat in the parameters ui. 

5 MDO mass distribution 

What do these new results tell us about the distribution of MDOs among galaxies? Let us assume initially 
that the MDO mass distribution of our sample depends only on x = M./Mb u i g c and is characterized by some 
other parameters u>; that is, that there is some function p(x | ui) dx which is the probability that a galaxy 
has an MDO with mass in the range [x, x + dx}. 

We experiment with several parameterizations P for p(x \ oj), as shown in Table 3. In each case, one of the 
parameters, /, is the fraction of galaxies with M. > 0, so that p(x | u>) is of the form 

p(x | u,P) = (1 - f)d(x) + fp+(x | uj,P), (11) 

where p+(x | lj,P) describes the distribution of MDOs with M, > 0. The N(lo) in p+(x | lj,P) is a 
normalizing factor chosen such that 



L 



P+ (x I u,P)dx = l. (12) 

o 

The parameterizations Ppl2 and Plg assume that there is a genuine ridge line in p(x) at x = xo, whereas 
the other three also test whether KR95's apparent ridge at x ~ 0.005 is just the upper envelope of some 
ridgeless p(x). 

For each parameterization P = (Ppli, PpL2> -Ps> -Pg> -Plg)) we nrs t see k the most likely set of parame- 
ters u> given our data D. By Baycs' theorem, the posterior distribution of u and mass-to-light ratios 
T = (Ti, • • • , T^) of the 32 galaxies is 

p(wT | D,P) oc p(w | P) p(T) p{D I ujY, P) 

(13) 



cx p(w | P) p(T) J p(D I Yx) p(x | to, P) dx, 



where ~p(D \ Yx) is a product of factors of the form of equation (8). The prior p(w | P) is assumed flat in 
the parameters u> given in Table 3. We are interested only in the parameters to, not in Y. Marginalizing the 
latter yields 

p(w | D, P) = J p(ujY I D, P) dT 



cx p{u | P) J p(D | x) p(x | u, P) dx, ^ 

N 



where we have defined 

p(D\x) = J p(D | Yx) p(T) dT, (15) 
7 



and similarly for p(Dj \ Xj). 



The posterior distributions p(w | D, P) for each parameterization are plotted on Figure 6. Table 4 lists the 
best-fitting parameters with their 68% confidence intervals and Figure 7 plots p(x | to, P) for the best-fitting 
parameters u> in each case. According to all parameterizations except Ppli, nearly all galaxies have MDOs 
(/ ~ 0.96) with means (x) ~ 0.01 and (logx) ~ —2.25, consistent with the KR95 interpretation. However, 
the best- fitting parameters from both PpL2 and Plg imply that there is a genuine ridge in p(x) at this mean 
x, whereas both Ppli and Ps say there is no ridge, since they prefer a < 0. The Gaussian parameterization 
Pq is inconclusive: there is not a strong lower limit on x in this case, since the most likely value of the 
other parameter A is comparable in size to x . 

Which of the five parameterizations gives the better description of the real p(x)? Using Bayes' theorem 
again, the plausibility of the parameterization P given the available data D is 



p(p|d) = p(p)pP|p) 



pig) 

P(D) 



J p(P | w,P)p(w | P) dw. 



(16) 



If we assume that all of the parameterizations are a priori equally likely, i.e., p(Ppli) = p(PpL2) = p(Ps) = 
p(P G ) = p(Plg), then we find that p(P PL2 | D) = 4.0p(P LG | D)~2x 10 5 p(P s | D) ~ 10 10 p(P PL i | D) ~ 
7 x 10 10 p(Pg I D): Ppl2 an d Plg provide by far the best description of the five. This result suggests that 
there really is a ridge in p(x) at logx ~ —2.2. 

It is also instructive to try to obtain a "non-parametric" estimate of p(x). We take n parameters uj\...u) n 
with n — 50. We define u>i as the probability that a randomly chosen galaxy has an MDO whose mass lies 
between Xi-i and Xi, where Xi runs logarithmically from x\ = 10~ 5 to xr M = 1, and xq — 0. A reasonable 
prior guess for p(x) (and therefore the w^) is a power law. So we choose 



, , /logw i+ i -21ogw i + logo; i _ 1 , 

iogP(-) = --TA ) , (17) 



i=2 



where the free parameter A controls how smooth (i.e., how far from a pure power law) we think an acceptable 
p(x) ought to be. We use 10 6 iterations of the Metropolis algorithm (Metropolis et al. 1953; see also Saha 
& Williams 1994) to obtain the posterior distribution p(w | D) for each uoi for a range of A. The results 
for A = 5 are plotted on Figure 7(b). The "non-parametric" distributions p(x) calculated in this way are 
broadly the same as the those obtained from the best parameterizations PpL2 and Plg- 

Thus far we have assumed that the MDO mass distribution depends only on x = M./Mb u i g o, but there is 
no good reason for assuming that M, should be correlated with the mass rather than, say, the luminosity of 
the bulge. Consider a more general form, 



Mq) \ L J V T , 

where a is a free parameter. The analysis above can be carried out with x replaced by x' . Setting a = 1 
tests the correlation of M, with Mbuige (the case we have just considered), whereas setting a = tests its 
correlation with L. The results of our calculations of p(P | D, a) for a range of a are plotted on Figure 8. 
Clearly, the case a = 1 is the most plausible - M. is much more strongly correlated with the mass of the 
bulge than the luminosity. 

Finally, a simple confirmation of the proportionality between M. and Mbuige can be obtained by splitting the 
sample of 32 galaxies into the most luminous half and the least luminous half, and then calculating p(aj | D, P) 
for each subsample for each parameterization. With the exception of the poorly fitting parameterization Ppli, 
we find that the best-fitting parameters u calculated using each subsample lie within the 95% confidence 
region of the parameters calculated using the full sample. 
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6 Conclusions 



We have examined a sample of 36 galaxy bulges and found that the kinematics of 32 of them are described 
well by two-integral axisymmetric models. Among these 32, a substantial MDO is required in all but four 
in order for our models to reproduce the observed kinematics. We have considered a range of models for 
the demography of these MDOs. In the best-fitting models about 96% of galaxies have an MDO. The mass 
of this MDO is strongly correlated with the bulge mass, with an MDO-to-bulge mass ratio of around 0.005. 
Possible explanations for this correlation have already been discussed by Faber et al. (1997). The galaxies 
without MDOs perhaps have a different formation history; one possible scenario has been put forward by 
Kormcndy (1993b). 

The mass-to-light ratios T fit by our models scale with luminosity L as T oc L - 2 , which is just the usual 
fundamental plane correlation. Since our models take full account of the shape of the light distribution of 
each galaxy, they rule out any attempts to explain the slope of the fundamental plane by a "non-homology" 
of the light profiles (e.g., Graham & Colless 1997, and references therein). Our models do not, however, 
consider the possibility of a systematic change in orbital anisotropy with luminosity (Ciotti et al. 1996). 

These results are based on an "assembly-line" approach to building galaxy models, which is necessarily less 
accurate than building models for each galaxy by hand. In particular: 

(i) Some or all of the galaxies may not be axisymmetric. 

(ii) Even if the galaxies are axisymmetric, our two-integral models are not the most general possible. For 
some or all of the galaxies, there may exist more general three-integral models that can reproduce 
the observed kinematics (and, indeed, the full line-of-sight velocity profiles) without needing to invoke 
MDOs - see Kormendy et al. (1997a) for an example. 

(iii) The selection criteria used to derive this sample are heterogeneous and impossible to quantify, although 
any biases introduced by properties such as luminosity, core size, and surface brightness are accounted 
for by the analysis procedure in §3. 

(iv) The assumption that the mass-to-light ratio is independent of position outside the centre may not be 
correct. 

Of the points above, the most important is perhaps (ii) - our conclusions are most uncertain due to our 
assumption of a two-integral distribution function. We should know soon whether more general three-integral 
models (e.g., Rix et al. 1997, Gebhardt et al. 1997, Richstone et al. 1997) will relax the need for MDOs in 
at least some of the galaxies in our sample. However, it is not yet clear what mechanism could effect just 
the right degree of radial anisotropy in each galaxy to cause the apparent correlation M, oc Mb u i ge as seen 
by our two-integral models. 
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A Appendix: Notes on individual galaxies 

M31: This galaxy has a double nucleus (e.g., Lauer et al. 1993), but axisymmetric models should still provide 
a reasonable description of the gross features of its kinematics. We use kinematical data from van der Marel 
ct al. (1994) and Kormcndy & Bender (1997). The former appear to measure major- axis radii from the 
photometric centre of the galaxy, rather than the kinematic centre, which we assume to be coincident with 
the fainter nucleus (e.g., Tremaine 1995). Thus we add 0.3 arcsec to van der Marel et al.'s quoted major-axis 
positions. We do not fit to kinematical data beyond 10 arcsec because the outer photometry we use (Kent 
1987) consists only of a major- and a minor-axis profile with no additional isophotc shape information. 
NGC 1600: Both Jcdrzejewski & Schechter (1989) and Bender, Saglia & Gerhard (1994) give major- and 
minor-axis kinematics for this galaxy. There are many more outlier points in the latter data, so we reject it. 
NGC 2778: Fisher, Illingworth & Franx (1995; FIF95) and Gonzalez (1993; G93) give major-axis pro- 
files. Both G93 and Jedrzejewski & Schechter (1989) give minor axis profiles. FIF95's central dispersion is 
inconsistent with the others, so we reject their data for this galaxy. 

NGC 3379: The best kinematical data comes from Gebhardt et al. (1997). We restrict our model fits to 
their ground-based data within 12 arcsec of the centre. 

NGC 4486: We use the blue G-band kinematics from van der Marel (1994), and reject his infrared kine- 
matics which are probably affected by template mismatch. In the same paper, van der Marel presents 
evidence that this galaxy is radially anisotropic in its outer parts. Therefore we restrict our model fits to 
the kinematics within the innermost 5 arcsec. 

NGC 4594: The ground-based outer photometry (Kormendy 1988) consists only of a major- and a minor- 
axis profile. Because of this, and because of the problems with dust obscuration, we only use kinematical 
data within 8 arcsec. 
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Figure 1. Major- and minor-axis projected second-moment profiles for a flattened (axis ratio q = 0.6) isotropic (rotating) Jaffc 
(1983) model viewed edge on. The curves show the (square root of the) classical second-order moments. The points plot the 
approximation (t> 2 + a 2 ) 1/2 where v and a are the parameters of the best-fitting Gaussians to the linc-of-sight velocity profiles. 
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Figure 2(a). Plots of \ 2 versus M, /M bulge for the four galaxies that our models do not describe well. The different curves on 
each plot correspond to different assumed inclinations. A reasonable fit would have \ 2 N dot ± (2A r dof ) 1 / 2 where the number 
of degrees of freedom, AT dof , is related to n, the number of kincmatical bins used in the fit, by N do f = n — 2. The heavy solid 
and dashed lines show x 2 = ^dof and x 2 = ^dof + (2./V do f) 1/2 respectively. 
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Figure 2(b). Kincmatical profiles of the best-fitting models along each slit position used for the galaxies in Figure 2(a). The 
plots show (v 2 + a 2 ) 1 ^ 2 (in units of kms -1 ) versus distance from the centre of the galaxy (in arcsec). The observed kinematics 
are plotted as circles, open or closed depending on which side of the galaxy the observation was made. The solid curves show 
the model predictions (convolved with the same seeing as the best ground-based observations) for the best-fitting value of M, 
and T for each galaxy. The results for each assumed inclination angle are plotted as separate curves. For comparison the dashed 
curves plot the model predictions with the same value of T as above but M, = 0. 
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Figure 3(a). The posterior distributions p(T, M 9 \ D) for all 32 galaxies that our models describe well. The vertical 
and horizontal axes are log(T/T Q ) and M,/M hu i sc respectively. Successive light contours indicate a factor of ten change in 
p(T,M. | D). The heavy contours enclose the 68% and 95% confidence regions on T and M.. 
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Figure 3(a). ..continued. 
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Figure 3(b). As for Figure 2(b), but for the galaxies that our models describe well. We also plot crosses to show the spatially 
binned, seeing-convolved model predictions for those cases where this quantity differs significantly from the unbinned model 
predictions (described by the curves). 
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Figure 3(b) ...continued. 
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Figure 3(b) ...continued. 
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Figure 4. The correlation between the MDO masses predicted by our models and those predicted by other methods. The 
crrorbars give the 68% confidence limits on M, . Sources for the other models are as follows: M31 — Richstone et al. (1990); 
M32 - van dcr Marel et al. (1997); NGC 3377 - Richstone et al. (1997); NGC 3115 - Kormendy et al. (1996a); NGC 4594 - 
Kormendy et al. (1996b); NGC 4486 - Harms et al. (1994). 
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Figure 5(a) and (b). The correlations between stellar mass-to-light ratio T and bulge luminosity L (left panel) and between 
MDO mass M. and M bu igc (right panel) produced by our models. The error bars give 68% confidence intervals. The solid lines 
plot Tm and M.,s t as described in the text (equations (9) and (10)). 
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Figure 5(c). The correlation of the residuals in the T-versus-L and M.-versus-Mb u ige fits. 
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Figure 6. The posterior distributions p(u> \ D, P) 
marginalized over / for (clockwise from above) P PL i, 
fpL2, fsj Pg and Plg- Successive light contours 
correspond to a factor of 10 change in p(oj | D, P). 
The heavy contours enclose the 68% and 95% con- 
fidence areas of the parameters oj. The most likely 
value of / within the 95% confidence area in all five 
cases is ~ 0.96. 



23 



-4 -2 
log x 







-4 -2 
log x 







Figure 7(a). (Left panel) The probability distributions p(x \ u),P) for the best-fitting parameters w. The heavy solid and 
dashed curves show results for Ppl2 and Plgi the two best-fitting cases. The lighter solid, dashed and dotted curves are for Ps, 
Ppli and Pq respectively, (b). (Right panel) The "non-parametric" probability distribution p(x) (heavy solid curve) and its 
68% confidence limits (heavy dashed curves) obtained using the Metropolis algorithm with A = 5. The rise in p(x) at small x 
is caused by those galaxies without an MDO. The best-fitting parameterized distributions Ppl2 and P l g are overlaid as the 
lighter solid and dashed curves respectively. 
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Figure 8. Variation of p(P | D,a) with a, where a is defined in equation (18). The solid and dashed curves show P = Ppl2 
and P = Plg respectively. The vertical scale does not extend down far enough to show the results for the other three 
paramcterizations . 
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Table 1. The galaxy sample. Column (1) gives the galaxy type: "S"=spiral bulge, "S0"=lenticular, "E"=elliptical; "P|"=cored, 
"\"=power law (Lauer et al. 1995). The absolute V magnitudes of the bulge or other hot component in column (2) are taken 
from Faber et al. (1997), and assume H a = 80 km s _1 Mpc -1 . i? m i„ in column (3) is the radius inside which we believe the 
galaxy light may be contaminated by non-stellar radiation. Column (4) gives the source of the outer photometry used (if 
available), otherwise it gives the effective radius R a u used for the outward extrapolation. Values of R„a obtained from the 
literature are followed by an '(1)', while an '(f)' follows those obtained by fitting to the HST photometry. Columns (5) to (9) 
list the kincmatical data used. For each exposure along each slit position, columns (5) and (6) give the position and FWHM 
of the seeing respectively. The maximum radius and the number of bins used by our models are given in columns (7) and (8). 
Finally, column (9) gives the source of the kinematical data. 
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Table 2. The best-fitting parameters (T, M, and x = M,/M bulsc ) with their 68% confidence intervals for the 32 galaxies that 
our models describe well. The assumed galaxy distance D and the luminosity L of the bulge or other hot stellar component are 
also listed. 
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Table 4. The best-fitting parameters ui and their 68% confidence limits for each assumed distribution p(x \ w, P). By definition 
< / < 1. The last two columns give the logarithm of the expectation value of x = M./M bu igo and the expectation value of 
log a; for those galaxies with M. ^ (both calculated from \> + (x | u>, P)). The mean (x) does not exist for Ppl2, while (log a;) 
docs not exist for Ppli- 
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